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Abstract The negative effective magnetic pressure instability discovered re- 
cently in direct numerical simulations (DNS) may play a crucial role in the 
formation of sunspots and active regions in the Sun and stars. This instability 
is caused by a negative contribution of turbulence to the effective mean Lorentz 
force (the sum of turbulent and non-turbulent contributions) and results in 
formation of large-scale inhomogeneous magnetic structures from initial uni- 
form magnetic field. Earlier investigations of this instability in DNS of stably 
stratified, externally forced, isothermal hydromagnetic turbulence in the regime 
of large plasma beta are now extended into the regime of larger scale separation 
ratios where the number of turbulent eddies in the computational domain is 
about 30. Strong spontaneous formation of large-scale magnetic structures is 
seen even without performing any spatial averaging. These structures encompass 
many turbulent eddies. The characteristic time of the instability is comparable 
to the turbulent diffusion time, L^/ryt, where 7]t is the turbulent diffusivity and L 
is the scale of the domain. DNS are used to confirm that the effective magnetic 
pressure does indeed become negative for magnetic field strengths below the 
equipartition field. The dependence of the effective magnetic pressure on the 
field strength is characterized by fit parameters that seem to show convergence 
for larger values of the magnetic Reynolds number. 
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1. Introduction 



The 1 1-year solar activity cycle manifests itself through the periodic variation of 
the sunspot number. S unspots consist of vertical magne tic fields with a strength 
of up to 3 kG (see, e.g. Parker . 1979t Ossendriiver , 2003 ). It is generally believed 
that these fields continue in a similarly concentrated fashion also beneath the 
surface in the form of magnetic flux tubes or fibers (,Priest . 1982.) . It was thought 
that fibral magnetic fields constitute a lower energy state and might therefore 
be preferred. Tube- like magneti c fields were f requentl y seen in hydroniagneti c 



turbulence simulations ((Nordlund et al. . 1992 : .Brandenburg et al. . 1995 , 1996I ). 



but those tubes were similar to the vortex tubes in hydrodynamic turbulence. 
Those vortex tubes have typical diameters comparable to the viscous length. 
Similarly, the aforementioned magnetic tubes have a thickness comparable to 
the resistive length. For the Sun, however, the resulting tube thickness would be 
too small to be relevant for sunspots. 

At the surface, strong magnetic flux concentrations also form in regions of 
strong flow convergence, but the size of these regions is too small for sunspots, 
because sunspots are usually much bigger than a single granular convection cell. 
In fact, a typical sunspot can have a diameter of some 30 pressure scale heights. 
The tremendous size of sunspots has therefore been used as an argument that 
they are not made near the surface, but at much greater depth far away from 
the surface. At the bottom of the convection zone the convection cells are big 
enough and could in principle be responsible for producing much bigger flux 
concentrations. Magnetic flux tub es can also form through the action of shea r 



as shown by various simulations (iCline et al.l . 120031 : iGuerrero fc Kapvlal . 120111) 



While shear is likely an important ingredient of the solar dynamo, it remains 
unclear whether the resulting magnetic tubes are really able to produce sunspots 
as a result of them piercing through the solar surface and, more importantly, 
whether one should consider them being tied to the deep shear layers of the Sun. 

In this paper we discuss an alternative scenario in which sunspots are shallow 
phenomena that are not anchored at the bottom of the convection zone. Various 
mechanisms have been discussed. Of particular interest here are mechanisms that 
are based on the sup pression of turbulence by magnetic fields. In the mechanism 
of ,Kitchatinov fc Ma zur (200^ it is the suppression of th e turbulent heat flux 
while in the mechanism of iRoeachevskii &: Kleeorin ( 2007 ) it is the suppression 
of the turbulent hydromagnetic pressure. Both mechanisms lead to a linear large- 
scale instability in a stratified medium. However, these mechanisms may be of 
different importance in different layers. 

In this study we focus mainly on the second mechanism, which has recently 
been studied in direct numerical simulations (DNS) as well as in mean-field 
calculations. This mechanism is called the negative effective magnetic pressure 
instability (or NEMPI for short). NEMPI is a co nvective type i nstability that is 
similar to the inter change instability in plasmas ( Tserkovnikovl 1960t Newcomb . 
Il96ll : IPriesd . 1982 ) and the magnetic buoyancy instabilitv (jParkeT 1966[ ). How- 
ever, the free energy in interchange and magnetic buoyancy instabilities is due 
to the gravitational field, while in NEMPI it is provided by the small-scale 
turbulence. NEMPI is caused by the suppression of turbulent hydromagnetic 
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pressure by the mean magnetic field. When fluid Reynolds numbers is larger than 
1 and the mean magnetic field is less than the equipartition magnetic fields, the 
negative turbulent contribution to the mean Lorentz force is enough large so that 
the effective mean magnetic pre ssure (the sum of tu r bulent and non-turbulent 
contr ibutions) become negative ( Kleeorin et al. . 19901 : Roeachevskii fc KleeorinL 



20071 ) ■ This is the main reason for the excitation of the large-scale instability 



that results in formation of large-scale inhomogcneous magnetic structures. 

The effect of the suppression of the turbulent heat flux has not yet been stud- 
i ed as extensivel y as NEMPI since the or i ginal paper of Kitchatinov fc Mazurl 
( 2000f ). although Kitchatinov fc Olemskov ( 20061) have used that model to study 
the decay of s unspots. Moreover, there is now some evidence that the effects 
anticipated bv lKitchatinov fc Mazur ( 2000l) may have already been in operation 
in various simulatio ns of solar convec tion where the spontaneous formation of 
pores ha s been seen (|Stein et aL , 2011 ) . Another example is the large-eddy simu- 
lation of Kitiashvili et al. ( 2010() where an imposed vertical large-scale magnetic 
field gets c oncen trated into giant vortices. This result is reminiscent of that of 



Tao et al.l (|1998[ ). who found a segregation into magnetized and nearly unmag- 
netized regions in stratified convection simulations. However, it has not yet been 
possible to obtain the large-scale magnetic structures resembling the sunspots, 
in models with a strong imposed mag n etic flu x tube structure at the bottom 



of the domain (see, e.g. Rempel et al. . 2009al lb: Cheung et al. . 2010t Rempel . 
2011allbl ). 



In the rest of this paper, we focus on NEMPI, which was first found in mean- 
field calculatio ns of a stra tified layer eeorin et al.l . ll996tlRogachevskii fc Kleeorinl [ 
[2007: Brandenburg et al.i . ,2010. ,201111 " However, th ose results remain ed uncon- 



vincing until NEMPI was also discovered in DNS (jBrandenburg et al.l . l2011a 



hereafter BKKMR). It it therefore most appropriate to begin our discussion 
with the latter. 



2. The model 

Following the earlier work of BKKMR, we solve the equations for the velocity 
the magnetic vector potential A, and the density p, 

= -c!Vp + JxB + p{f + g) + V- {2vpS), (1) 

^ = UxB + rjV'A, (2) 

^ = -V-pU, (3) 

where u is the kinematic viscosity, 77 is the magnetic diffusivity due to Spitzer 
conductivity of the plasma, B = Bq + V x A is the magnetic field, Bq = 
(0,-Bo,0) is the imposed uniform field, J = V x B/po the current density, 
Po is the vacuum permeability, = \ + C^j.i) ~ ^^ij'^ ' ^ is the traceless 
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rate of strain tensor, and commas denote partial differentiation. The forcing 
function / consists of random, white-in-time, plane non-polarized waves with a 
certain average wavenumber. The forcing strength is such that the turbulent rms 
velocity is approximately independent of z with Umis = (it^)^^^ ~ 0.1 Cg. The 
gravitational acceleration g ~ (0,0, —g) is chosen such that kiHp = 1, which 
leads to a density contrast between bottom and top of exp(27r) « 535. Here, 
Hp = Cg/g is the density scale height. We consider a domain of size Lx ^ Ly X 
in Cartesian coordinates (x, y, z), with periodic boundary conditions in the x and 
y directions and stress-free perfectly conducting boundaries at top and bottom 
(z = -iiLz/2). The volume-averaged density is therefore constant in time and 
equal to its initial value, po = {p)- 

Our simulations are characterized by the scale separation ratio, fct /fci , the fluid 
Reynolds number Re = Ui-ms/i/fcf, the magnetic Prandtl number Pr^f = v/rj 
and the magnetic Reynold s number Rej\f = RePr7\/. Following earlier work 



([Brandenburg et al.l l2011bf ) , it is clear that NEMPI is more effective for small 
values of PrM, so here we choose Prj\f = 0.5 and Rcm in the range 0.7-70. 
The magnetic field is expressed in units of the local equipartition field strength 
near the top, — y^JIopuTms^ 

while Bq is specified in units of the volume 
averaged value, Boqo = ^/MoPo^rms- In addition to visualizations of the actual 
magnetic field, wc also monitor ABy = By — Bq, where By is an average over y 
and a certain time interval At. Time is expressed in eddy turnover times, Tto ~ 
{urinsk[)~^ ■ For comparison, we also consider the turbulent-diffusive timescale, 
Ttd — ivtokf)"^, where rjto = Wims/3fcf is the estimated turbulent magnetic 
diffusivity. Another diagnostic quantity is the rms magnetic field in the Fourier 
mode of = fci, referred to as Bi, which is here taken as an average over 
2 < kiz < 3. and is close to the top at kiz ~ tt. (Note that Bi docs not include 
the imposed field Bq at k = 0.) 

The simulations are performed with the Pencil CodeQ which uses sixth- 
order explicit finite differences in space and a third-order accurate time stepping 
method. We use numerical resolutions of 128'^ and 256'^ mesh points when Lx = 
Ly = Lz; and 1024 x 128^ when Lx — 8Ly = 8Lz- To capture mean-field 
effects on the slower turbulent-diffusive timescale, which is rtd/Tto ~ 3kf/kf 
times smaller than the dynamical timescale, we perform simulations for several 
thousand turnover times. 



3. Results 

In simulations, the clearest indication of a spontaneous flux concentration is seen 
when the scale separation ratio is large. In BKKMR, we only used k{/ki = 15. 
Here we consider calculations where this ratio is twice as large. A useful diag- 
nostic is the magnetic field averaged along the direction of the imposed field, 
i.e., along the y direction. In particular, wc shall be looking at the y component 
of the field, i.e., we look at {ABy)y/ B^q- To sec the effect even more clearly, we 



^http : //pencil-code ■ googlecode . com | 
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Figure 1. Visualizations of By{x, z,t) for different times. Time is indicated in turbulent 
diffusive times, {T]tokf)~'^ , corresponding to about 5000 turnover times, i.e., t = SOOO/nrmsfcf . 
RcM = 18 and Pim = 0.5. 

perform an additional time average over about one hundred turnover times. This 
average is then referred to as {ABy)yt. In Figure [T] we show {ABy)yt/ Bcq for 
k{/ki = 30. The other parameters are Rcm = 18 and PrM — 0.5. An inhomoge- 
neous magnetic structure forms first near the surface (at t/rtd ~ 0.79), but then 
the structure propagates downward. This is consistent with our interpretation 
that this is caused by negative effective magnetic pressure operating on the 
scale of many turbulent eddies. Indeed, a local decrease of the effective magnetic 
pressure must be compensated by an increase in gas pressure, which implies 
higher density, so the structure becomes heavier and sinks in the nonlinear 
stage of NEMPI. This is also seen in three-dimensional visualizations without 
averaging; see Figure [2] with the same parameters as in Figure [TJ 

To confirm that NEMPI really is a large-scale instability, we would expect 
to see an exponential growth phase. This is shown in the right-hand panel of 
Figure [3J where we show the growth of Bi versus time. We give time both in 
turbulent-diffusive times (lower abscissa) as well as in eddy turnover times (upper 
abscissa) . We do indeed see that there is an exponential growth phase which lasts 
for about one turbulent-diffusive time, i.e., the growth rate is comparable to 
(,Vtok'^)~^ , where rjt « Urms/S/cf is the expected turbulent magnetic diffusivity. 
However, compared with the eddy turnover time, (urms^f)"^, the turbulent- 
diffusive time scale is 3(fcf/fci)^ times slower. This illustrates that NEMPI is 
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Figure 2. Visualizations of By on the periphery of the domain for different times, indicated 
in turbulent-diffusive times, for the same run as in Figure [T] 




Figure 3. Time dependence of -Bmis (left panel) and Bi (right panel) for the same run as in 
Figure[T] where fcf/fci = 30 and Bo/^cqO = 0.05. 

indeed a very slow process compared with, for example, the saturation of the 
overall rms magnetic field (left-hand panel of Figure [3]) . 

While the rate of NEMPI may be too slow to explain the relatively rapid 
appearance of sunspots on a time scale of days, it would explain how one can 
produce magnetic structures on length and time scales much bigger than the 
naturally occurring scales in the upper layers of the Sun. This discrepancy is 
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Figure 4. Time dependence of i?rma (left panel) and Bi (right panel) for a run used in 
BKKMR with fcf/fci = 15 and Bo/BeqO = 0.05. 



exactly one of the reasons why one normaUy places the formation of active 



regions at the bottom of the convection zone (IGolub et al.l . Il98l[ ). Here we see 



a clear example that this conclusion may be not correct. 

Next, we turn to a simulation where the scale separation ratio, kf/ki, is only 
half as big, i.e., kf/ki = 15. In that case we also see an exponential growth phase, 
but the growth is slower (even in terms of turbulent-diffusive times), lasts longer, 
and amplifies Bi only by about one order of magnitude; sec Figure SI 

Again, in agreement with related work involving the suppression of turbulent 
heat flux, the most unstable mode has a horizontal wavelength comparable to 
the vertical scale height of the layer; see Fig. 2 of iKitchatinov &: Mazur 
hereafter referred to as KM. This is also the case for NEMPI, as is shown in 
FigureO where we compare an instantaneous plot of {ABy)y{x, z) with a time 
averaged one, {ABy)yt, making the appearance of large-scale structures even 
more pronounced. Next, in Figure |6] we compare cross-sections of ABy{x) (for 
fixed values of j/, z, and t; top panel) with corresponding y averages (middle 
panel) and yt averages (bottom panel). Here, Rcm = 36, fcf/fci = 15, and 
Ba/Bcqo = 0.05. Without averaging, no clear magnetic structure is seen yet, 
but the structures become clearly more pronounced with y and t averaging. 
Runs with similar parameters have been shown in BKKMR, for a computational 
domain whose x extent is 2-k /ki instead of 167r/fci. 

The large-scale flux concentrations have an amplitude of only Bi « 0.1 Boq 
and are therefore not easily seen in single snapshots, where the field reaches 
peak strengths comparable to i?oq- Furthermore, as for any linear instability, 
the flux concentrations form a repetitive pattern, and are in that sense similar 
to flux concentrations seen in the calculations of KM that were based on the 
magnetic suppression of the turbulent heat flux. However, there are indications 
that at larger values of Rej\/ , flux concentrations occur more rarely, which might 
be more realistic in view of astrophysical applications. 
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Figure 5. Visualization of By{x,z) for an elongated box (1024 X 128^ mesh points) with 
Kc]\{ = 36 at a time during the statistically steady state. The top panel shows the y average 
{ABy)y/Bcq at one time while the lower panel shows an additional time average {ABy)yt / Bcq 
covering about 80 turnover times. 
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Figure 6. x dependence of the field for an elongated box (1024 X 128^ mesh points) with 
Rcjv/ = 36 at kiz = 2 showing ABy/Beq at j/ = (top panel), its y average {ABy)y / Bcq 
(middle panel), as well as an additional time average {ABy)yt/ Beq (bottom panel) covering 
about 80 turnover times. The dash-dotted line gives the level of the imposed field. 
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4. Quantifying the negative effective magnetic pressure effect 

An impo rtant c onditio n for the formatio n of structures by tlie meclianisms of 



KM and iRoeacIievskii fc KleeorinI (j2007f ) is sufficient scale scale separation. In 



other words, both the suppression of convective heat flux and the suppression 
of total effective pressure have in common that they only work if a substantial 
number of eddies is involved in a turbulent structure under consideration. If this 
is not the case, turbulent transport coefficients become increasingly inefficient. 
This is a natural feature of large-scale effects such as this one. Here, we mean 
by turbulent transport coefficients any of the mean-field coefficients that relate 
correlation functions of small-scale quantities in terms of large-scale quantities. 

Small- and large-scale quantities refer here simply to a suitably defined aver- 
aging procedure so that the velocity f7, for example, can be split into a mean 
(or large-scale) quantity U and a fluctuation u = U U. An important example 
of a turbulent transport coefficient is the turbulent viscosity that emerges when 
relating the Reynolds stress Uiuj to the spatial derivatives of the mean flow U 
in the averaged momentum equation, 

Q _ Q 

Here, p is the average density, and correlations with density fluctuations are 
neglected. The simplest parameterization for UiUj is 

y-u] = -vtiUi^j + U ■j.j) - i-h5tjV ■ U, (5) 

where Vt is the turbulent viscosity and fit is the turbulent bulk viscosity. This 
relation is also known as the Boussinesq ansatz, especially when contrasted with 
representations where the A effec t is included, w hich is responsible for producing 
differential rotation in the Sun (Rudiger, 19891 ). Here, however, we shall focus 



on magnetic effects. 

In general, when there are magnetic fields, the right-hand side of equation of 
motion Q has to be replaced by the sum of Reynolds and Maxwell stresses 

H,^ = puTirj - hbj/fiQ + ^b'^/fio, (6) 

where the superscript f indicates contributions from the fluctuating field. Fur- 
thermore, in the presence of a mean magnetic field, symmetry arguments allow 
one to write down additional components, in particular those proportional to 

SijB^ and BiBj. 

— f 

Expressing H^^ in terms of the mean field, the leading terms are 

= qsB,Bj/fin - igpJyB^/^o + (7) 

where the dots indicate the presence of additional terms such as the turbulent 
viscosity term mentioned earlier and terms that enter when the stratification 
affect the anisotropy of the turbulence further. Note in particular the definition 
of the signs of the terms involving the functions qs{B) and qp(B). This becomes 
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clear when writing down the mean Maxwell stress resulting from both mean and 
fluctuating fields, i.e., 



-{l-qs)BiBj/nQ + ^{l-qp)S,jB^/flO 



(8) 



Thus, the signs arc defined such that for positive Qs and qp the effects of magnetic 
stress and magnetic pressure are reduced and the signs of the net effects may even 
change. E quations Q and jH) have been derived using the spect r al r r elaxation 



approach (jKleeorin et all . ll990L Il996l: iRogachevskii fc Klc eorinl l2007l ) and the 



renormalization procedure (jKleeorin fc Rogachevskii . ,19941 



A broad range of different DN S in stratified turbu lence (jBrandenburg et al 



2011al lb ) or turbulent convection ( Kapvla et al. , 2011 ) have now confirmed that 
Qp is positive for Re a/ > 1, but Qs is small and negative. A positive value of qs 
(but with large error bars ) was originally reported for unstratified turbulence 
([Brandenburg et al.Ll2010l). _jater , stratifi ed simulations with isothermal stable 
strat ification (Brandenburg et al. . l2011bl ) and convectively unstable stratifica- 
tion ( Kapvla et al.l . 2011 ) show that it is small and negative. Nevertheless, qp{B) 
is consistently larger than unity provided Re a/ > 1 and above unity while B/ B^q 
is below a certain critical value that is around 0.5. This is shown in Figure [7l 
where we plot the effective magnetic pressure, 



7'eff(/3) = ^[1- gp (/?)]/?' versus I3 = \B\/B 



cq 



(9) 



for different values of Rcm using PrM = 0.5 and k[/ki = 15. Note that the 
minimum of Ves{P) is deeper for the case with Rej\f = 11 and becomes then 
shallower. 

It should be noted that the coefficients for the fit formulae arc chosen such 
that the plot is most accurate for small values of /?, because this is what is most 
important for the onset of NEMPI. Especially for larger values of Rejvf the fit 
formula does no longer reproduce correctly the critical value /? — /3crit for which 
VcR changes sign. 

Note that /3crit is well below unity. This implies that it is probably not possible 
to produce flux concentrations stronger than half the equipartition field strength. 
So, making sunspots with this mechanism alone is maybe unlikely, and other 
effects such as that of KM may be needed. It is possible that this mechanisms 
works preferentially in the uppermost layers, provided enough flux has already 
been accumulated. This may then be achieved with NEMPI which also works in 
somewhat lower layers. 

To compare the resulting functio ns VesiP) in a system atic fashion for different 



parameters, we use the fit formula (jKemel et al.l . I2012i 



0i 



1 



/3? 



where 



(10) 



The resulting dependencies /3p(ReM), /3*(ReM), and (/po(ReM) are shown in 
Figure [HI We see that /3^(ReM) varies relatively little between 0.1 and 0.2 and is 
typically around 0.15. For small values of Rcj\/, /3p(ReA/) drops from 1 to 0.1 and 
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Figure 7. Normalized effective magnetic pressure, "Poffl/S), for low (upper panel) and higher 
(lower panel) values of Re^vj . The solid lines represent the fits to the data shown as dotted 
lines. 



stays then approximately constant, while gpo(ReM) rises proportional to Rcj^^ 
for Rca/ < 10 and then levels off at a value around 40. 

The significance of a positive qs value comes from mean-field simulations with 
Qs > indicati ng the formation of three-dimensional (non-axisymmetric) flux 



concentrations ( Brandenburg et al. . 2010l ). This result was la ter identified to 



be a direct consequence of having gs > (jKemel et al Before making 

any further conclusions, it is important to assess the effect of other terms that 
have been neglected. Two of them are related to the vertical stratification, i.e. 
additional terms in Eq. ([7]) that are proportional to giQj and QiB-j + QiBj with t r 



being gravity. The coefficient of the former term seems to be small (jKapvla et al 



20111 ) and the second only has an effect when there is a vertical or inclined 



imposed magnetic field. However, there could be other terms such as JiJj as 
well as JiBj and JjBi when the scale separation is not enough large. 
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Figure 8. Re a/ dependence of /3p, /3*, and gpo for Prm = 0.5 and kf/k\ = 15. 



5. Conclusions 



In the present paper we have performed detailed investigations of NEMPI de- 
tected recently by BKKMR. Most notably, we have extended the values of the 
scale separation ratio, k[/ki, from 15 to 30. In this case, the spontaneous for- 
mation of magnetic structures becomes particularly evident and can clearly be 
noticed even without any averaging. Whether or not the particular structures 
seen in DNS really have a correspondence to phenomena in the Sun, cannot be 
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answered at the moment, because our model is still quite unrealistic in many 
respects. For example in the Sun, fcf and Urms change with depth, which is not 
currently taken into account in DNS. Also, of course, the stratification is not 
isothermal, but cony ectively unstable. However, DNS in turbulent convection by 



Kapvla et al.l (j201lh have shown that Vcs{l3) still has a negative minimum in 
that case. 

With regards to the production of sunspots, it is likely that NEMPI will shut 
off before the magnetic energy density has reached values comparable with the 
internal energy of the gas, as is the case in sunspots. Thus, some other mechanism 
is still needed to push the field of fiux concentrations into that regime. One likely 
candidate is the mechanism of KM where the suppression of convective heat fiux 
by the magne tic field is crucial . This impression is further justified by resent 



calculations of IStein et al.l (|201lh where pores are seen to form spontaneously in 
a simulation where horizontal magnetic fields are injected at the bottom of the 
domain. 

Pores are small sunspots, so something else is needed to make these struc- 
tures bigger and to amplify this mechanism further. Again, the answer could be 
related to scale separation. Thus, it will now be important to undertake detailed 
investigations of instabilities in strongly stratified layers with finite heat fiux and 
finite magnetic field. 
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